A FIRST PRINCIPLES STUDY OF CU CLUSTERS 
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Abstract 

In this communication we study the equihbrium shapes and energetics of Cu clus- 
ters of various sizes upto 20 atoms using the Full-Potential Tight Binding Muffin-tin 
Orbitals Molecular Dynamics. We compare our results with earlier works by physicists 
and chemists using different methodologies. 

1 Introduction 

Many of the first principles molecular dynamics approaches to the study of clus- 
ters depend upon the construction of suitable pseudopotentials for the constituent 
atoms. Transition metal clusters require perhaps, alternative treatments ( p8| , ^ 
^). The deep potentials associated with their d-orbitals are not particularly 
amenable to the pseudopotential approach. In this communication we shall describe 
a study of Cu clusters using a full-potential LMTO based molecular dynamics. 
Experimentation on the electronic and cohesive properties of transition metal 

^ m m m m c 



clusters have been extensive jg, ^, |2S|, 0, |3|, |5T|, §g, 0, 0. The 

smaller Cu clusters have been exhaustively studied by quantum chemists ^, ^ 

H, m, H m, m, 1^, O. An excellent early 



io|, y, yj, y, yj |i7|. 



review of the field has been made by Ozin [34|. The main issues addressed in these 



works were : whether small clusters had characteristics of bulk metals and in what 
way they differed from them in respect to cohesive energies, ionization potentials 
and magnetism. Recently Apai et al ^ conducted EXAFS studies of Cu clusters 
supported on carbon. Similar studies of Au and Ag clusters were carried out by 
Balerna et al [|16l and Montano et al |E9i. These studies indicate, as one would 



expect for these metals, that the localized d- electrons play an important role in the 
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electronic structure. Hence, the d-states and their interactions with the extended 
s-states need to be carefully accounted for in a proper theoretical treatment of these 
materials. 

We may classify the theoretical approaches into five groups : 

(i) In the first group are the Hartree-Fock and Xa descriptions ([§§!). In this class 
we have the the self-consistent-field-Xa (@), the ab-initio self-consistent- 
field (SCF) using model potentials (|]l5l) and the SCF with relaxation effects 

M). 



(ii) In the next group are the methods based on the local spin density (LSDA) 
both without and including self-interaction corrections (|T^). Salahub and 
coworkers have argued that it is essential to include the gradient corrections 
in the LSDA in order to treat clusters properly ([^, ^) since 
the bonding charge density in a small cluster is highly inhomogeneous. 

(iii) The third group includes tight-binding type methods. These include the ex- 
tended Hiickel methods (0, |12|, [1^), re-parameterized Hiickel with the Wolfsberg 



Helmholtz approximation for the off-diagonal terms (IQ) ^^"^ those with more 
flexible forms for them ( ) . This group also has the linear combination of 
atomic orbitals based SCF methods (pSf). We also have the empirical tight- 
binding (TB) or the Linear Combination of Atomic Orbitals (LCAO) methods 
(p7|, p8|). These are at best quahtative, since the assumption of transferabil- 



ity of the Hamiltonian parameters is definitely of questionable validity. 

(iv) In the fourth group we have the effective potential methods which include 
the embedded atom pair and many-body potentials (|^^ |6^) and the effective 
medium theory ([0, |4^) with one-electron correlation included ([1^, |2^, |78 



The equivalent crystal theory (ECT) (|j72|, ^ |7^) also belongs to this class of 



empirical potential methods and is capable of dealing with very large clusters. 

(v) Finally we have attempts at using the tight-binding linearized muffin-tin or- 
bitals (TB-LMTO) method coupled to simulated annealing. In the appli- 
cation of this method to clusters there are several outstanding problems. The 
treatment of the interstitial region outside the muffin-tin spheres centered at 
the atomic positions is difficult. Unlike the bulk, where the interstitial region 
is small and inflating the muffin-tins to slightly overlapping atomic spheres can 
do away with the interstitial altogether, for clusters this is certainly not so. 
As the atoms move about, the atomic spheres may not overlap and the inter- 
stitial contribution is signiflcant. One may try to overcome this by enclosing 
the cluster with layers of empty spheres carrying charge but not atoms. This 
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complicates the actual calculations and the justification of extrapolation of the 
TB-LMTO parameters beyond the 5% range on either side of the equilibrium 
value is not valid. 

A number of molecular dynamics studies of the geometrical and electronic struc- 
ture of small clusters of various elements ([^2|, ^) have been performed. The 
ab-initio molecular dynamics (MD) approach developed by Car and Parinello [21 
(CP) has been one of the most promising developments in this area. The method is 
based on the pseudopotential technique and therefore faces problems when dealing 
with the rather localized Delectrons of transition metals. Efficient soft pseudopo- 
tentials for transition metals are still not available and the CP generally is never 
applied to transition metal clusters. 

Simple alkali metals clusters are fairly well described by the spherical jellium 
model. The quasi-free valence electrons occupy single-particle states in an effective 
spherically symmetric box potential. This is rather insensitive to the geometry of 
the atomic arrangement inside the cluster. Consequently one obtains a pronounced 



shell closing effect ([0). Although the noble metals Cu, Ag and Au have closed 
d-shells and singly occupied outermost s-shell structures and several authors have 
suggested that there should be a close similarity to the shell closing effect in simple 
alkali metals, cohesive studies in the bulk metal and a series of EXAFS studies of 
Cu clusters supported on carbon (p, [Q) indicate that the d-electrons through 
their hybridization with the s-electrons play an important role in the electronic 
structure and binding energy of these systems. have also indicated through a 
series of experiments which include mass spectroscopy, oxygen and water absorp- 
tion, that there is a competition of jellium-like electronic behaviour and icosahedral 
geometrical closure effects in small copper clusters. 

In this chapter, we shall turn to the molecular version of the full-potential lin- 
earized muffin-tin orbital two-centre-fit (TCP) method suggested by and 



to carry out an ab initio study of Cu clusters ranging in size between 10 and 20. 

2 The two-center fit method : TCF 

The molecular version of the full potential-LMTO two-centre-fit (TCF) method 
utilizes the philosophy of Muffin-Tin Orbitals methods. It is based on the Den- 
sity Functional Theory in the Local Density Approximation. The electron-electron 
interaction is treated approximately. In practice : 

+ Veff{r)]i;,{r) = eitPi{r) (1) 
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where, 



and, 



p{l) = l^/il^id) 



Veff = VN{r) + 2jj^^^dY + l^.c{p{r)) 

The first step is the solution of the Schrodinger equation in a very unpleasant 
potential with Coulomb singularities. As in most approaches we use the variational 
approach. We choose a basis of representation {0m(ll)} such that 

m 

The problem reduces to a matrix eigenvalue problem : 

He — eSc 

Computation effort scales as ~ (matrix dimension)^. Our approach tries to use 
a minimal basis set at the expense of a rather complicated formulation. The basis 
is built up of Hankel functions Hj^ diverging at r = R^^ augmented inside the 
muffin-tin spheres by solutions u(r)Y2,(f) of the Schrodinger equation : 

(^ + 1) 



<{r) = 



+ y(r) 



UL{r) 



with boundary conditions such that its logarithmic derivative matches that of the 
Hankel function. Any matrix element in this basis then can be written as : 



'^iL\0\(t)jL>) 



I + I 

JSk JI 



(t^:L{r)0(t>,v{r)dh 



(2) 



The Hankel functions associated with a muffin-tin at can be written in terms 
of a Bessel function at Rj as YjL" SiL,kL"JkL" the structure matrix S depends entirely 
on the geometric arrangement of the muffin-tins. The first integral becomes : 



^^^^SiL^kL"SjL',kL"{JL"\0\JL")sk + X! ^3L',iL{Hi\0\JL)si- ■ 

+ E SiL,iL'{JL\0\HL)s, + {Hl\0\Hl)s, 

k=j,yti 



(3) 
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These are easy to calculate and there is a separation of atomic and structural infor- 
mation. 

Most of the interstitial integral can be obtained from the muffin-tin spheres by 
using the fact that, in the interstitial, the basis are solutions of the Helmholtz 
equation, and using the Green theorem : 

Jl Ki K-2 1^ JSk 

(PI (- V^) = kI cl^lhdh (4) 

If the potential here is a constant we can get bye with the above. But for clusters 
this is definitely not so. In the molecular FPLMTO we use a tabulation technique. 
We expand the product : 

<i>Kr)h{r) = Y^Ciixmim) 

m 

where XmiL) is another set of muffin-tin centered Hankel functions. In practice we 
put two atoms along the z-axis and make accurate numerical expansion by least 
squares fit for different distances and tabulate C^(d) : 

j^XmiL)Xn{r)d^r 
l(t>*{r)4>j{r.)Xm{r)d^r 

This is the two-centre fit table (TCF). For arbitrary geometry then we may easily 
calculate the necessary matrix elements by a fitting procedure to the table. The 
procedure id fast. 

For molecular dynamics, the problem arises from the fact that the Pulay terms 
in the force are impossibly difficult to calculate directly as the basis set changes in 
a complicated manner when atoms move. To do the molecular dynamics, we use 
the Harris functional procedure as follows : At a time step tq we obtain the self- 
consistent charge density p(r, tq) using the FP-LMTO procedure. At a neighbouring 
time To+T we hazard a guess Pg{r, tq + r) and obtain 

E(r) = EH[pg(r,ro + T)] 

= E [Veff] - J Pg{r)V,ff [P.(r)] + U [p,(r)] + E,, [p,(r)] 



A 
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To find the force on an atom, we simply move tlie atom witli its surrounding 
cliarge density in a given direction. Tlie force is given by 



dE 

For dynamics we use the Verlet algorithm 



m 



where n denotes the time step of length At. We can now either do straightforward 
molecular dynamics, but this often leads to unphysical heating/cooling of the system 
if our time steps are too large. For small time steps the procedure is inordinately 
slow. We add an extra friction term carefully F_ ^ K — "ymr. Methfessel and 
Schilfgaarde |Q have also used a free dynamics with feedback to overcome the 
above difficulty. 



3 Results 

We have chosen the various parameters for the FP-LMTO based on optimizing 
results for bulk Cu and the dimer. The values of were chosen from optimum 
bulk calculations. The muffin tin radii were chosen as 1.9 Ato produce the bond 
length and binding energies of the Cu dimer correctly. For augmentation within the 
sphere we have used 4s, 4p, 3d, 4/ and 5g functions (^maa;=3). For representation 
of interstitial functions we have used five values with angular momentum cutoffs 
i^ax = 4,4 6,2 and 1. 

The optimum bond length was determined by varying the dimer bond length 
from 4.1 to 4.2 atomic units and calculating the total energy at each bond length. 
We found the optimum bond length to be 4.16 a.u with the binding energy (B.E) 
equal to 1.469 eV/atom. The table ^ lists the various theoretical and experimental 
values for the bond length and binding energies per atom. It is well known that while 
the Hartree-Fock tends to under-bind, the LDA over-binds. Our bond lengths should 
then be smaller and binding energies larger than experimental values. This is borne 
out by the table. Clearly both the self interaction correction (SIC) and the gradient 
correction (GGA) improves matters. The TB-LMTO value of 0.23 eV/atom (||52||) 
is much too low and probably indicates serious lacuna in the treatment of clusters 
in that work rather than in the TB-LMTO itself. 
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Table 1: Bond lengths and Binding energies/atom for Cu2 dimer 



Bond length 


Binding Energy 


Method 


Reference 


(a.u.) 


(eV/atom) 







Theoretical Results 



4.55 




ab initio SCF 


N 




4.28 


0.923 


LSD Pseudo. 






4.20 


1.025 


LSD Pseudo. + SIC 


[7S 




4.28 


0.975 


LSD Pseudo. + CI 


[ZS 




4.56 


0.34 


Hartree-Fock 


i 




4.55 


1.04 


CI 


i 




4.20 


1.50 


X-a 


i 




4.17 


1.30 


LCGO-DFT 


IS] 


4.30 


1.13 


LCGO-GGA 


m 


4.16 


1.369 


FP-LMTO-TCF 


Our work 




0.23 


TB-LMTO 


m 




Experimental Results 


4.195 


0.99 


Expt 


m 




4.21 


1.03 
1.04 




i 
m 
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Table 2: Ionization potentials for Cu2 by various methods 



IP (eV) 



Method 



Reference 



5.65 ab initio SCF 

6.04 ab initio SCF Koopman 

7.987 LSD Pseudopotential 

8.237 LSD Pseudopotential + SIC 

6.37 Hartree-Fock 

7.37 Modified Hartree-Fock-I 

7.89 Modified Hartree-Fock-II 

7.35 X-a 

7.64 EHT 

5.70 LCAO-SCF 

8.69 LCGO-DFT 

7.904 LCGO-DFT + CCA 

8.22 FP-LMTO-TCF 



m 
m 
m 



i 
i 
m 

m 

Our work 



The table ^ shows the ionization potential (IP) for Cu2 dimers, calculated as the 
difference between the total energies of neutral Cu2 and the ion, using various 
methods. The experimental values quoted range between 7.904±0.04 quoted by 
Calamici et al and 7.37 of Joyes and Leleyer It is quite clear that for 

the smaller clusters the generalized gradient corrections (GGA) to the local density 
approximation is very important ([|19|)- Our FP-LMTO does not incorporate the 
GGA and hence leads to slightly larger values of the IP. The importance of self- 
interaction corrections (SIC) is not clear for dimers. Wang includes SIC and 
obtains a higher value of the IP. Our work does not include the SIC. 

The first test of the predictability of various methods first appear for Cu-^. The 
accompanying figure 1 shows the lowest energy structures predicted for the trimer. 
Miyoshi et al |5^ find both the structures (O) and (A) to be almost degenerate in 
energy. The vertex angles are found to be 77.6° for (O) and 51.7° for (A). Calamici 
et al \T^\ find the structure (O) to be most stable with vertex angle 66.86° without 
SIC and 66.58° with SIC. The other structure (A) lies 0.023 eV higher in energy. 
Wang |79| finds the obtuse triangle shown on the right to be the stable structure. 
This has a vertex angle of 162°. He concludes that the SIC correction is essential and 
finds the acute triangle with a vertex angle of 47° to be the most stable. However, 
even with the SIC the structure quoted is rather different from other methods. Our 
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Table 3: Bond lengths, binding energies and ionization potentials of Cus 



Side length 


Binding Energy 


IP 


Reference 


(a.u.) 


(eV/atom) 


(eV) 





Theoretical Results 



4.72 


1.018 


4.39 


4.50 


0.942 


7.144 


5.41 


0.753 


6.018 


4.35 




6.33 


4.28 


1.34 


6.46 


4.45 


1.12 


5.795 




0.68 




4.30 


1.598 


6.40 



m 
m 



m 
m 



52 



Our work 



Experimental Results 



1.02 



5.48±0.5 



prediction agrees reasonably well with the structure (O) of Miyoshi et al and 
(O) of Calamici et al [|l^. The vertex angle is 65° in our case. The isosceles shape 
is expected because of the possible Jahn- Teller distortion in Cu^,. 

The table |^ compares the bond lengths, binding energies and ionization energies 
of he Cu^ trimer. We find the binding energy per atom to be 1.598 eV/atom which 
is higher than that for the linear configuration by 0.124 eV/atom. Over-binding 
because of the LDA is again observed. The ionization energy drops for the trimer 
and regains its value again for CM4. This has been observed in all the earlier works 
quoted and in experiment. 

For N=4 we find the rhombus starting structure to lead to the most stable 
structure followed by the square and the tetrahedron in decreasing order of stability. 
Our prediction matches exactly with that of Akeby et al P[ and Calamici et al [P!U 



who also predicted the sequence rhombus, square and tetrahedron. The larger 
rhombus angle turns out to be 120° which agrees well with the prediction of 122° 
by Calamici et al [0. Our ionization potential is 7.90 eV, which agrees not badly 
with 7.0±0.6 eV found experimentally. The TB-LMTO ([^) predicts the order of 
stability to be the tetrahedron, the rhombus and the square in decreasing stability. 
This does not match with any other work and possibly has its origin in the problem 
talked about earlier. 

For Ciis we find the trigonal bipyramid with B.E. 2.187 eV/atom to be the 
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most stable structure followed by the square pyramid where the difference in B.E. 
between the two structures is .056 eV only. Akeby et al also obtain the trigonal 
bipyramid to be more stable than the square pyramid agreeing with our calculations. 
Calamici et al finds another structure, the flat pentagonal trapezoid to be almost 
degenerate; actually 0.009 eV lower in energy than the trigonal bipyramid. They 
find the square pyramid to be more than 0.309 eV higher in energy. We would like 
to emphasize with Calamici et al that for the smaller clusters the GGA may play a 
crucial role in stabilizing certain structures. 

Figure 3 shows the variation of the ionization potential with cluster size. The 
troughs at n=3 and n=5 agree well with earlier works as well as experimental results 

O). 

For N=6 we have considered two starting structures the square bipyramid (octa- 
hedron) and the capped trigonal bipyramid which is obtained by capping one face 
of the trigonal bipyramid so that the capping atom is equidistant from all the three 
atoms on the face. We find the capped trigonal bipyramid to be the most stable 
structure with bond energy equal to 2.405 eV/atom which is 0.040 eV/atom higher 
than the square bipyramid (octahedron). The TB-LMTO calculations predict the 
octahedron to be the most stable structure compared to other random structures. 
Also the numerical value of 1.56 eV/atom for the octahedron obtained from the 
TB-LMTO calculations is much lower compared to our value. 

The pentagonal bipyramid, the capped square bipyramid and the bicapped trig- 
onal bipyramid were considered as the starting structures for our calculations for 
N=7. We find the pentagonal bipyramid to be the most stable structure in ac- 
cordance with Akeby et al . but at variance with the TB-LMTO results. The 
bicapped trigonal bipyramid is slightly higher in energy (0.002 eV/atom) than the 
capped square pyramid in our calculations. 

For Cuq we considered three starting structures as shown in the table of which 
the capped pentagonal bipyramid turns out to be the most stable followed by the 
bicapped square bipyramid and the cube. TB-LMTO predicts the antiprism fol- 
lowed by the bi-tent structure and the cube. Both the methods find the cube to 
b ethe least stable though our B.E. for the cube is 0.592 eV/atom higher than the 
TB-LMTO results. 

In the case of Cug, we considered the tricapped square bipyramid and the bi- 
capped pentagonal bipyramid with the capping atoms on adjacent and non- adjacent 
faces. The tricapped square bipyramid was found to be the most stable structure 
followed by the bicapped pentagonal bipyramid with the capping atoms on adja- 
cent faces (lower by only 0.006 eV/atom) and the bicapped pentagonal bipyramid 
(non-adjacent faces) lower by 0.025 eV/atom than the most stable structure in this 
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range. The stable shapes for 6 < N < 9 are shown in figures 4. 

Figure 5 shows a plot of the binding energy versus cluster size for N=2 to 9. 
The relative stability of the clusters {2E{N) - E{N + 1) - ^^(A^ - 1)) is also plotted 
on the same graph. Cluster sizes N=3,5 and 8 show up as more stable. This has 



been predicted experimentally earlier by Knickelbein Katakuse |^ has also 
observed N=8 to be a stable structure in their experimental observations. 

Figure 6 shows a plot for the HOMO-LUMO gap versus cluster size for the most 
stable clusters. In the theoretical results of Akeby et al |0] the HOMO-LUMO gap 



for the Cu% cluster is determined to be 1.93 eV while hammers and Borstal 
report a value of 1.91 eV. Our calculated value for the HOMO-LUMO gap for the 
most stable structure (capped pentagonal bipyramid) is 1.156 which is lower than 
both the reported values. The HOMO-LUMO gap does show a peak at N=8 in our 
calculations but we cannot conclude from this point that this is a manifestation of 
shell closure. We also see a minimum in the HOMO-LUMO gap value at N=6 unlike 
in . Moreover pronounced odd-even alterations in the HOMO-LUMO gap values 
as predicted by the shell model (||75|) are not recognizable in our calculations. 

The N=10 cluster shape is a close competition between the tetracapped trigonal 
bipyramid which is obtained by capping the N=9 cluster on another face and the 
structure shown on the right hand side of figure 6. Our calculations indicate that 
the former is more stable, however, the energy difference is smaller than the errors 
involved in the FP-LMTO itself. From N=ll to N=13 the clusters grow towards 
the stable icosahedron. These shapes indicate that probably our prediction is valid. 

For = 12 we started from a configuration which is an icosahedron with a void 
at the centre. Rapidly the structure evolved to the icosahedron with one exterior 
atom removed. 

For = 13 we studied carefully two possible structures : the cubo-octahedron 
(shown on the left in figure 9 and the icosahedron, shown on the right of the same 
figure. Our calculations indicate that even if we begin with the cubo-octahedron as 
our starting structures, the cluster rapidly settles down to the icosahedron. Earlier 
Valkealahti and Manninen |]78| had also used effective medium-molecular dynamics 
and shown that the cubo-octahedron is unstable and rapidly changes over to the 



stable icosahedron. Winter et al [pil] have argued from experimental observations 



that the shell structure seen in the smaller clusters is overshadowed by icosahedral 
closures from A^ = 13 onwards. 

For A^ = 15 and A^ = 16 we see near-degenerate structures. The lower-energy 
structure has atoms on neighboring faces of the icosahedron. There is also another 
structure, differing in energy by about 1%, in which the "extra" atoms are on 
non-neighboring faces of the icosahedron. For A^ = 17 the two different starting 



11 



structures both anneal to an icosahedron with four atoms on neighboring faces. 

The = 19 has a very stable structure : the double icosahedron, further con- 
firming the conjecture of Winter et al |^ regarding icosahedral closure. For = 20 
The equatorial addition was found to be more stable by about 1%. We expect as 
the size increases, the cluster structure becomes more spherical. Note that we see 
no evidence for the very open structure reported to have been obtained by hammers 
and Borstal ||5^ for = 20 through simulated annealing. 

Figure 12 shows the binding energy and the homo-lumo gap for the clusters 
iV = 11 to = 20. We note that the signatures of shell closure we observed in 
the smaller clusters is overtaken by geometric closures and the icosahedron based 
closed structures are the more stable. 
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FIGURE CAPTIONS 



1 Various shape predictions for the Cus trimer 

2 Stable configurations for CU4 and Cus 

3 Variation of the ionization potential with cluster size 

4 Stable configurations for Cug to Cug 

5 The binding energy per atom and its curvature for Cu2 to Cug 

6 The homo-lumo gap for Cu2 to Cug 

7 The structures for Cuio 

8 The stable structures for Cun and Cuu 

9 The structures for Cuia 

10 The structures for Cui4 - Cuig 

11 The structures for Cuig - CU20 

12 The binding energy per atom and the homo-lumo gap for Cun to Cu2o 
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